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1 INTRODUCTION 



Eviden ce for the existence of dark matter dates back to IZwickvl 
dl933h with studies of the kinematics of galaxies in the Coma 
Cluster, which required the presence of a massive, smooth and 
dark component generating the cluster gravitational potential. More 
recently, astrophysical and cosmological observations as well as 
simulations and theoretical arguments have provided further in- 
dication for the existence of dark matter, and hopes that it may 
be detected directly in particle accelerators (see for example 
iFrandsen et al and references therein). These developments 

include the flatness of galaxy rotation curves (e.g., Bosmal [T978I ; 



iBosma & van der Kruitffl97a ; iRubin et alJll978l . TT980n . the mass 



of galaxy cluster s inferred either by th eir X-ray emission (e.g., 
lAllenetalJl2002l:IVikhlinin etal.|[2006h or by gravitational lens - 
ing (e. g. lClowe et aT1l2006l : lUmetsu et al.|[201ll : Mira et alj|201 ll) . 
the acoustic oscillations measur ed in the cosmic mi crowave back- 
ground (see WMAP papers, e. g Jarosik et al. 20 iT j) and in galaxy 
surveys (e. g. lEisenstein et al . 2005; Sanchez et al. 20125, and de- 
tailed studies of structure format ion on numerical sim u lations with 
ever increasing precision (e.g. ISpringel et al.l 120051; lAlimi et all 
I20TI . 

In the standard scenario, dark matter is assumed to be com- 
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ABSTRACT 

We use the stacked gravitatio nal lensing mass prof ile of four high-mass (M> 1O 15 M ) galaxy 
clusters around z ~ 0.3 from lUmetsu et all (1201 ll) to fit density profiles of phenomenological 
(NFW, Einasto, Sersic, Stadel and Hernquist) and theoretical (non-singular Isothermal Sphere, 
DARKexp and Kang & He) models of the dark matter distribution. We account for large- 
scale structure effects, including a 2-halo term in the analysis. We find that the Stadel model 
provides the best fit to the data as measured by the reduced % 2 . It is followed by the generalized 
NFW profile with a free inner slope and by the Einasto profile. The NF W model prov i des th e 
best fit if we neglect the 2-halo term, in agreement with results from lUmetsu et all d201 lb . 
Among the theoretical profiles, the DARKexp model with a single form parameter has the 
best performance, almost identical to that of the Stadel profile. This may indicate a connection 
between this theoretical model and the phenomenology of dark matter halos, shedding light 
on the dynamical basis of empirical profiles which emerge from numerical simulations. 

Key words: cosmology: dark matter; galaxies: clusters; galaxies: halos 



posed of particles that interact gravitationally but not electromag- 
netically. Within this picture, simulations of structure formation 
have shown a number of interesting results regarding the final states 
of systems of gravitationally interacting particles. For example, it 
has been shown that dark matter is cold, i.e. its particles must have 
had non-relativistic velocities around the epoch of recombination, 
otherwise structures like galaxies and galaxy clusters would be 
more diffuse than they appear due to free-streaming. 

Another feature that emerges from simulations, and is con- 
firmed by observations, is that these systems seem to achieve a fi- 
nal st ate of equilibrium, displayin g nearly universal density profiles 
p (r) l lNavarrp et al.lll99"ol Il997l) and pseudo phase-space profiles 
dl2C 



p (r) /a 3 (r) jTavlor & Navarroll200lh where a (r) is either the ra- 
dial or total velocity dispersion. This is intriguing because these 
non-collisional systems interact only through gravitational effects, 
making them quite different from e.g. a molecular gas in a box. We 
are then led to approach the question of how a non-collisional pro- 
cess could bring a gravitational system to an equilibrium state in a 
time scale of the age of the Universe. 

In an attempt to answer this kind of question, iLvnden-Belll 

i 19671) developed the mechanism of violent relaxation, in which 
the system's constituents, e.g. stars or dark matter particles, interact 
mainly with a time-varying average gravitational field, for which 
the time scale to achieve an equilibrium state is many orders of 
magnitude smaller than that of two-particle interactions. However 
there are issues in this approach, like infinite masses and mass seg- 
regation. These happen because the model generates density pro- 
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files similar to that of an Isothermal Sphe re (see below ) . For m ore 
details on thes e and other critics t o the iLvnden-Belll dl967t) ap- 
proach, see e.g. lHiorth & Williams! 11201 On and references therein. 

Regardless of the discussions about the statistical process in- 
volved, the density profiles provided by simulations are observed in 
real objects like galaxy clusters, in which the role of the baryonic 
component is relatively small compared to that of dark matter. Phe- 
nomenologically, one can argue that if simulations provide density 
profiles which match those in observed data, that means that the 
assumptions made in the simulations are likely correct, and the ob- 
served features are consequences of the gravitational interaction. 
Nonetheless, a deeper understanding of the physical mechanisms 
that lead to equilibrium in gravitational systems is certainly desir- 
able. In fact, ignoring this issue would be equivalent to making sim- 
ulations of a molecular gas in a box and computing gas pressure and 
velocity distribution from the simulated results, with no regards to 
the kinetic gas theory developed by Maxwell, Boltzmann and oth- 
ers. 

In order to have a better dynamical picture of gravitational sys- 
tems in general, and of dark matter halos in particular, first principle 
models have been developed to explain the features seen in simu- 
lations and observations. In particular, there has been a great effort 
to make predictions of the three-dimensional density profile p (r) 
of dark matter halos. The connection with observations is made via 
the surface density profile E (R) projected in the line-of-sight jcm 



T.(R) = J <£t|| p{x\\,R) 



(1) 



where r = (xm ,R) and R is the projected distance on the plane 
of the sky. Compared to galaxies, for which dissipative effects of 
cold baryons are important, galaxy clusters are excellent to test the 
distribution of dark matter, because in clusters most of the baryons 
are hot and dissipate less. Thus, the total density profile, inferred 
e.g. from gravitational lensing measurements, provides reliable in- 
formation about the dark matter density profile. In fact, lensing is 
particularly interesting in the determination of the observed density 
profile of galaxy clusters, because it does not require assumptions 
of hydrostatic equilibrium, as in dynamical methods. 

In this paper we use the stacked surface density profile from 
four massive galaxy clusters with similar mass and redshift to test 
both phenomenological and theoretical models for their density 
profiles. We only consider spherically symmetrical models. In ^2] 
we briefly describe the cluster data used in this work. In §|3]and $4] 
we present the phenomenological and theoretically motivated mod- 
els tested. In ^5]we present the halo model, which allows us to in- 
clude large-scale structure effects on the observed profiles. Our re- 
sults are presented in ^6] and discussed in ^7] When necessary, we 
use the values of WMAP 7-year data release d-Farosik et alj|201 if) 
for the cosmological parameters. 



2 DATA 



We use the data of lUmetsu et al.l d201 ll) . who combined weak- 
lensing shear, magnification, and strong-lensing measurements of 
four high-mass (M > 1O I5 M ) galaxy clusters (A1689, A1703, 
A370, C10024+17) with redshifts around z « 0.3. The strong lens- 
ing data was based on Hubble Space Telescope observations for 
the central regions of those clusters (typically, R < 150 kpc//z), 
and combined with independent weak-lensing data obtained by 
lUmetsu et al.l J201 lb . extending to the outer regions (R < 3.5 
Mpc//7) of the clusters. The surface density profiles of the four clus- 
ters were stacked in order to reduce the cosmic noise and smooth 



effects due to asphericity or presence of substructures. Hereafter 
we assume that the radial shape of the mean density profile ob- 
tained i n this way is rep resentative of dark matter halos in equi- 
librium dGao et al]l2012l) . For m ore observational inf ormation on 
individual clusters, see Table 1 of lUmetsu et"aHd201ll) . 



3 PHENOMENOLOGICAL MODELS 

A number of phenomenological models for density profiles of 
dark matter halos and galaxy clusters have been proposed as 
parametrized functions that fit reasonably well simulations and ob- 
servations, with no regards to fundamental principles or theoretical 
motivation. Below we make a brief description of the models that 
we test. 



3.1 NFW profile 

The NFW profile was proposed by iNavarro et all (1 19961 1 19971) in 
order to fit the data of N-body cold dark matter (CDM) simulations, 
after stacking many halos. It is given by 



PW 



(r/r s )(l + r/r s y 



(2) 



where p., and r s are scale parameters. It often represents the best fit 
model to observed data of galaxy clusters; this would also be the 
case in this work if we did not include large-scale structure effects 
in the analysis, as discussed below. The NFW profile has an ana- 
lytica l expression for the surface density profile E (R) itBartelmannl 
1 19961) . given by 



T.(R)=2p s r s F(R/r s ), 



(3) 



where 



F(X): 



1 



X 2 -l 
1 

3' 

1 



2 ll-X' 
-arctanh \ , 

l+X 



VT^xi 



l- 



2 IX- 1 

=arctan i 



Vx^T 



X+l 



(X<1) 

(X = l) 
(X>1) 



The NFW profile has a non-physical divergence at the origin, vary- 
ing as in the inner regions. In its outer parts it varies as 
implying another unrealistic property of an infinite total mass. One 
way to circumvent the latter divergence is to truncate the profile at 

a maximum radius, e.g. the virial radius. 

A common g eneralization of the NFW profile dZhaol fl99ol 
Ijing & Sutoll2000l) is obtained by setting the inner slope as a free 
parameter a (for NFW, a = — 1): 



(4) 



{r/r s T{\ + r/r s f- al 

Following lUmetsu et al.l ( feoill) we will refer to this generalized 
model as gNFW. 

3.2 Sersic profile 

The Sersic profile (see ISersidl 19631) was proposed in order to fit 
the light distribution in spheroidal galaxies and has been also used 
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to fit simulated data iMerritt et al.l2006h . It is defined as a projected 
surface density profile that has the form 



£ (R) = L e exp { -b n [(R/R e ) l/ " - l] } 



(5) 



where Z e is the surface brightness at the effective radius R e and b„ 
is a function of n obtained by imposing that the luminosity inside 
the effective radius is half the total luminosity. The r elation betw een 
b n and n is well approximated by b n ~ In — 0.324 lCiottill99lh . 



3.3 Einasto profile 

The Eina sto profile is a t hree-dimensional version of the Sersic pro- 
file (see lEinastdfl965» . It was proposed to describe the surface 
brightness of elliptical galaxies. Recently it has also been used to 
fit data from N-body CDM simul ations, giving results comparable 
to the NFW profile in some cases ( Navarro et al. 2004; Merritt et al.l 
l2005tl2006l : lGao et alJl2008l ; lNavarro etai1l2010l) . It is given by 



p(r) = p. s exp{-2,7 [(r/r^) 1 /" - l] } , (6) 

where p. s and r_2 are scale parameters. lLapi & Cavalierel fcoioh 
discuss a possible dynamical basis that could generate a profile 
for which Einasto is a good approximation. In their Appendix A, 
Mam on et al. have obtained a polynomial approximation to 

better than 0.8% for the expression of the surface density, in the 
intervals 3.5 < n < 6.5 and —2 < log 10 {R/r^i) ^ 2. 

3.4 Stadel profile 

This profile was proposed to fit sim ulated data of galaxy-size dark 
matter halos (see IStad el et al. 2009)). It has the form 



p(r)=poexp|-X[ln(l+r/r s )] 2 | , 



(7) 



which resembles somewhat the Einasto profile and similarly gives 
a finite density po at the origin. It can also be written as 

Po 



P(r) 



\\n(l+r/r s ) : 



(8) 



and in this way it resembles power-law profiles. Noticing that the 
shape param eter assumed almos t the same value X = 0. 1 in different 
simulations. [Stadel et alj J2009T) proposed to fix this parameter and 
promote the model into a two-parameter profile. Here we let A. be a 
free parameter and obtain a different value for it. 



3.5 Hernquist profile 

The lHernguist1 l ll990h profile has the functional form 



(9) 



(r/r s )(l + r/r s y 

and differs from the NFW profile only in the outer parts, where it 
varies as r~ 4 . It was proposed, not as a fit to simulated or observed 
data, but because it provides analytical expressions for dynamical 
quantities, such as the gravitational potential, the energy distribu- 
tion function, the density of states as well as the surface density, 
which is given by 



Z(R)=2p s r s G(R/r s 



where 



G(X) 



[(2 + X 2 )H pQ-3] 
2(1-X 2 ) 2 



(10) 



(11) 



and 



H(X): 



1 



i, 

i 

— — < 

Vx^T 



sech _1 X 



(X<1) 
(X = l) 
(X>1) 



which implies that lim G(X) = 2/15. 

After addition of the 2-halo term explained in ^5] the phe- 
nomenological profiles described above are shown in F i g. [2] a long 
with the galaxy cluster stacked data from lUmetsu 



4 THEORETICAL MODELS 

Some of the theoretical models we investigate here are based on 
the hypothesis of hydrostatic equilibrium between the gravitational 
attraction and the pressure P(r) due to velocity dispersion in an 
isotropic distribution: 



dP 
~dr 



-P{r 



GM(r) 



(12) 



where p (r) is the mass density profile and M (r) is the total mass 
inside radius r: 



M(r) = J' dr'4nr l2 p (/) . 
Combining Eqs. d!2t and J 1 3b we have 

-4ltGr 2 p (r) . 



d 

7/7 



r 1 dP 



p(r) dr 



(13) 



(14) 



Choosing the equation of state P(p) determines the model, 
and Eq. d 14b can then be (numerically) solved to give the density 
profile p (r). 



4.1 (non-singular) Isothermal Sphere 

The (non- singular) Isothermal Sphere is based on the equation of 
state of an ideal gas P = nk^T , which locally becomes 



P{r) = p(r), 

m 



(15) 



where m is the mass of the constituent particle. Using Eq. ( 1151 in 
Eq. (O we have 



d 2 p / dp\ 2 dp 



dr 7 



dr 



2p—+4nGXrp 3 =0, 
dr 



(16) 



where X = m /kgT . This represents a particular case of the so-called 
Lane-Emden equation. The (non-singular) Isothermal Sphere has 
null slope at the origin and oscillates around the Singular Isother - 
mal Sphere (p oc r~ 2 ) for large radii feinnev & Tremainell2008l) . 
Thus we solve Eq. ( 116l > numerically imposing the boundary condi- 
tions dp/dr(0) = and p (0) = po, where po is a free parameter. 



4.2 Kang & He models 

The entropy per unit mass s r of an ideal gas, written as a function 
of pressure p r and density profile p (r) is 



r = m( P rp- 5/2 ), 



(17) 
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and the Jeans equation is written as 



Pr 



dr r 



GM(r) 



(18) 



where (3=1— ^o" 2 , + o" 2 , j / (2<7 2 ) is the velocity anisotropy param- 
eter, written in t erms of the velocity dispersions in the three spher- 
ical coordinates. lKang&Hd l f2011bh define a generalized pressure 
Pas 



dr dr r 
and a phenomenological entropy as 

, = ln(V 3 / 2 p- 5/2 



(19) 



(20) 



such that the resulting system of equations is independent of p\ This 
effectively reduces Eq. d!8) to Eq. d!2t . Using the variational prin- 
ciple, the entropy per unit mass, Eq. ( 1201 . is then used to maximize 
the total entropy 5 



Anr 2 psdr = J 47tr 2 pln 



3/2.-5/2 



dr. 



(21) 



subject to the constraints of conservation of total energy and the 
virial theorem. This procedure results in the following equation of 
state 



p = lP + /jP y , 



(22) 



where y = 3/5. We will refer to this model as "Kang & He". The 
constant X is a Lagrange multiplier and /j is an integration constant, 
both related to total mass and energy of the system. This equa- 
tion of state reduces to that of an ideal gas Eq. H5) f or fj = and 
X = m/ ksT. Following a similar but different approach. Ka ng & Ffej 
obtain the same equation, but now with y = 4/5. We will 
call this last model "Kang & He 2". 

In order to use the equation of state Eq. d22t in Eq. d!4t . we 
need to solve for P(p) in Eq. d22t . so as to turn Eq. d!4t into an 
equation for p(r). This is n ot possible for gener al values of y, so 
after differentiating Eq. d22l . iKang &~Hel ( 1201 lbh propose approx- 
imating P ~ p/K obtaining 

dP _ 1 dp 

Ir ~ 



(23) 



X + y/j(X/pY- y dr ' 

Inserting Eq. d23t into Eq. d!4t , one obtains a second order differ- 
ential equation for p (r), which can be numerically solved imposing 
again dp/dr(0) = and p(0) = p . 

It is possible to follow a different approach, inserting p(r) 
from Eq. d22t into Eq. d!4t . thus obtaining an equation for P(r). 
After solving this equation numerically for P(r), p(r) can be ob- 
tained from Eq. d22t . This approach proves to give slightly better 
results (although similar to the original Kang & He's ) in the fitting 
procedure, so that is what we used. 

4.3 DARKexp 

The DARKexp model faiorth & Williams! l20ld : 

IWilliams & Hiorfhl 1201 d) is significantly different from the 
previous models, because it does not take into account a possi- 
ble equation of state to be used in the hydrostatic equilibrium, 
Eq. d 14b - Instead, it deals with statistical mechanical arguments 
to (indirectly) derive the distribution function and determine the 
density profile. 

For the discussion below, let us define a dimensionless density 
p = p/po and a dimensionless distance x = r/a, where po and a are 



scale parameters. The particle's energy per unit mass E = <!> + v 2 /2, 
where $ is the gravitational potential and v the particle velocity, can 
be written as 



9- 



2v 2' 



(24) 



where v ? = \J a 2 PoG and we defined the positive and dimensionless 
quantities e = —E/v g 2 and (p = — <J>/v ? 2 . 

The DARKexp model is based on two main assumptions. 
First, because dark matter in halos is collisionless, it is argued that, 
after the system reaches an equilibrium, each particle retains its in- 
dividual energy, and thus a Boltzmann-like function must be used, 
not in the distribution function /(e) (average number of particles 
per state of energy e), but in the number of particles per unit en- 
ergy N (e) oc f (e) g (e), where g (e) is t he den sity of states (number 
of states per unit energy); see iBinnevI dl982l) . The other feature of 
the model is that it properly considers the possibility of low occu- 
pation n umbers, which results in a cutoff similar to that of King 
models i lKing|ll96d : lMadseiJ[l99l) . These two features imply that 
the number of particles per unit energy e must be given by 



N(e) = eV°- e -l, 



(25) 



where (po is the shape parameter representing the central potential. 

In models that predict the distribution function /(e), the den- 
sity profile is obtained after int egrating over all possible velocities 



sity profile is obtained alter inte 
(see iBinnev & Tremaine1l2008t) 

47t J 



h, v 2 J^-/(e) 



(26) 



where / (e) is considered to be dimensionless. The equation above 
also assumes that the velocities are isotropic, contrary to what is 
seen in simulated ACDM halos (e. g. Lemze et alj2012h and in ob- 
servational analysis (e. g. iBiviano & Katgertll2004h . With the help 
of Eq. d24b . we have 



p (x) = An / def (e) y/2[<f(x)-E] . (27) 
Jo 

The density profile is finally obtained by solving Poisson's equation 

-47tp (x) . (28) 



V 2 cp(x) 



However, if the model predicts N(e), as in the case of the 
DARKexp, we need to use an iterative approach. Here we follow 
the procedure o f lBinnevlfl982h . We start by guessing an initial esti- 
mate of the density profile p (x) and calculate the resulting potential 

as 



<p (jt) = 47t - f dx'x' 2 p (x') + [°° dx'x'p (Y) 

_ X JO J x 

Next, we compute the density of states as 

g(E) = {4lt) 2 J drr 2 J dvv 2 8 Qv 2 + 
which in terms of the dimensionless quantities results in 



g{e) = 16% 



a v g J dxxr y 2 [<p (x) — e] , 



(29) 



(30) 



(3D 



where x max is such that (p (x max ) = e. We then use the N (e) defined 
in the model, Eq. d25l >. to compute the dimensionless distribution 
function 



/(e)=« 3 v g 



g(e) 



(32) 
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Finally, we use Eq. J27t to obtain a new p (jc) and iterate the pro- 
cess. We find that after about 20 iterations the model converges to 
a density profile independent of the initial guess. 

Fig. [3] shows the theoretical models described above, after 
adding to them the 2-halo ter m explained in jf5l al ong with the 
galaxy cluster stacked data from lUmetsu et al.ldioill) . 



5 HALO MODEL 

When considering cluster profiles that extend to sufficiently large 
radii, large-scale corrections must be taken into account. For dark 
matter halos of a given mass M and redshift z, the halo-mass corre- 
lation function, defined as %hm{ r ) = {S/i( x )°m( x + r ))> represents 
the excess density of matter at a distance r = |r| from the halo 
center, i.e. it is a measure of the average observed halo profile 
(PobsM): 



(PobsW) 



(33) 



The halo model (see lCoorav & Shetbll2002l) allows us to esti- 
mate cosmological correlations from the properties of dark matter 
halos, seen as the building blocks of cosmic structure. In this con- 
text the halo- mass correlation function is given by a sum of two 
contributions dHavashi & Whitel2008USchmidt et alj|2009h 



%hm{r) 



P\h{r) 



-^(M)^(r). 



(34) 



Here Pi/,(r) represents the 1-halo contribution or true halo 
profile from matter within the halo itself; this is the term described 
by all models presented in ^3]and ^4] The second term on the right- 
hand side of Eq. d34b represents the 2-halo contribution from the 
large-scale structure of the Universe, given by the linear matter cor- 
relation function fjjj(r) and the linear halo bias b^(M). 

Projected lensing measurements are sensitive to the average 
observed overdensity 8p [, s (r) = {p bs( r )) — pm- Therefore, com- 
bining Eqs. l |33l > and l |34t we find 



5PobsM = P\h(r) + p2h(r), 
where the 2-halo term is given by 
P2f,(r) = p m b^(M)% L m (r). 



(35) 



(36) 



The observed surface density profile Z b s (7?) at projected dis- 
tance R is obtained using Eq. |QJ: 



EobsW = J dx l} Sp obs (x h R)=I, lh (R)+Iah(R) 



(37) 



where Y.u,(R) is defined from Pi/,(r), and simi larly for T.2i,(R). We 
estimate bu(M) from the fit to simulations of iTinker et al.l J2OI0I) 
and %m( r ) as me Fourier transform of the linear matte r power spec- 
trum pL(k) obtained from CAMB jLewis et alj2000l) . 



2;t 2 / 



dkk 2 pt(k)- 



kr 



(38) 



Finally we assume a flat Universe with cosmological param- 
eters set by the b est fit values of WMAP 7-year data release 
jjarosik et alfeoi ll) . 

In Fig. [T] we illustrate the effect of the 2-halo term for mas- 
sive halos. We assume a 1-halo profile given by the NFW model 
from 33.11 and consider halos with average virial mass Af VI > = 
1.56 x 10 15 M Q //? and redshift z = 0.32. For this redshift and cos- 
mology, th e virial overdensity relat ive to the mean matter density 
is A ~ 263 tervan & Normanlll998h , and the Tinker fitting formula 



NFW 

2-halo term 

NFW + 2-halo term 




R(Mpc/h) 



Figure 1. NFW profile (dotted), 2-halo term correction (dashed) and the 
sum of both (solid) representing the final best fit are shown, along with data 
points from lUmetsu et al] J201 ll) . 



gives bi(M263,z) = 10.98. The NFW model underestimates the ob- 
served profile for r > 1 Mpc//z, where the 2-halo term becomes in- 
creasingly important. We add the computed £2/1 to the 1-halo mod- 
els of 5j3]and ^4]before fitting them to data. Since the 2-halo contri- 
bution depends only on the fixed cosmology, this does not introduce 
any extra parameter. 



6 RESULTS 

We fit the various mo dels to the data with t he help of the Minuit 
package developed by Ijames & Rood Jl975l) . We compute the % 2 , 
i.e. the minimum value of Q, given by 

(39) 



where 



A, =E r (Ri)-L D (Ri) 



(40) 



Y.T{Ri) is the surface density from Eq. |[T) for a given model 
evaluated at radius /?,, (/?,■) is the surface mass density from 
lUmetsu et alj d201ll) lensing analysis of the data, and Vij is the er- 
ror covariance matrix between data points i and j. The data consist 
of 15 correlated points. 

Hereafter, the metric we use to compare the various models is 
the % 2 per degree of freedom, or reduced x 2 , defined as 

,2 

(41) 



r 

V 



where v = 15 — N p is the number of degrees of freedom given 15 
data points and N p parameters. 

We note that a more rigorous statistical analysis would involve 
the correct calculation of the integral of the % 2 distribution. Such 
detailed analysis is beyond the scope of this work, since we are 
mainly interested in investigating what class of models provide a 
reasonable description of the data, as inferred from a simple rank- 
ing criterion. We believe though that our conclusions would remain 
unchanged if a more detailed statistical analy sis was employed. 

Fig.|2]shows the data points obtained bv lUmetsu et alj fcoill) 
and the best fits for all the phenomenological models discussed 
above, after the addition of the 2-halo term represented by the or- 
ange dashed line in Fig.[T] 
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Figure 2. Best fit cur ves and %v f° r the ph enomenological models along 
with data points from lUmetsu et all bOl ll) . The bottom panel shows the 
relative difference between the best fits of each model and the Stadel profile. 
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Figure 3. B est fit curves and %v for the theoretical models, along with data 
points from lUmetsu et al.l 120111) . The bottom panel shows the relative dif- 
ference between the best fits of each model and DARKexp. 



Profile 




Xv 


Shape parameter 


Stadel 


3 


0.465 


?i = 0.25 ±0.04 


gNFW 


3 


0.471 


a = 0.74 ±0.44 
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Hernquist 
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NFW 
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0.484 




Sersic 


3 


0.504 
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Table 1. Fit results for the phenomenological models. The column N p in- 
dicates the total number of model parameters, Xv shows the reduced y} de- 
fined in Eq. 14 U and the last column shows the best estimate for the shape 
parameter of the model. 

The Stadel profile, with 3 parameters and Xv = 0.465, rep- 
resents the best fit, with best value of the shape parameter X = 
0.25 ± 0.04. It is followed by the generalized version of NFW, 
gNFW (3 parameters), with %v = °-471 and a = 0.74 ± 0.44. 
Next is the Einasto profile (3 parameters), with Xv = 0.475 and 
n = 3.80±0.61. Then we have the Hernquist model (2 scale param- 
eters) with %y = 0.482, followed by the NFW model (2 parameters) 
with Xv = 0.484. Finally, for the Sersic profile we obtain x v = 0.504 
and n = 2.42 ± 0.34. These results are summarized in TableQ] 

Fig. [3] shows the fits for the theoretical models. The best fit 
was obtained for the DARKexp model with x v = 0.468. In order to 
generate this model, we did 25 iterations of the procedure described 
in § 14.31 with 10 5 logarithmic bins in r. The best fit value for the 
shape parameter was (po = 3.00±0.48. The other three models, the 
Isothermal Sphere and the 2 variants predicted by Kang & He, were 
generated in 10 6 logarithmic bins in r. 

Among these last three models, the best fit is for the "Kang & 
He" model, with y = 3/5, for which x v = 2.350 with X = Xc 2 = 

(5.44±0.10) x 10 4 and ( u = J u(c 3 /po) 2/5 = 13.85±0.03 for the 



Profile 




Xv 




Shape parameter 


DARKexp 


3 


0.468 




(p = 3.00±0.48 


KH 


4 


2.350 


I 


= (5.44±0.10) x 10 4 
fi= 13.85 ±0.03 


Isoth. Sph. 


3 


2.603 


X 


= (5.76±0.16) x 10 4 


KH2 


4 


2.670 


I 


= (5.29±0.38)x 10 4 
£ = 346 ±268 



Table 2. Fit results for the theoretical models. Columns defined as in Ta- 
bled 

shape parameters, where c is the speed of light in vacuum and po 
is the scale parameter for the density. The Isothermal Sphere gives 
Xv = 2.603 and X = (5.76 ±0.16) x 10 4 for the shape parameter, 
followed by "Kang & He 2", with y = 4/5, for which Xv = 2.670 

with X= (5.29 ±0.38) x 10 4 and fi = /j{^/Po) 1/5 = 346 ±268. 
These results are summarized in Table [2] 

6.1 Neglecting the 2-halo term 

We have also considered the results of fitting the models without 
adding the 2-halo term. These fits are summarized in Figs.|4]and|5] 
and Tables[3]and|4]for the phenomenological and theoretical mod- 
els respectively. 

In this scenario, the NFW profile produces the overall best 
fit, with Xv = 0.449 and the gNFW profile results in Xv = 0.474, 
with the best fit val ue a = 0.89 ± 0.37. T hese values are identical 
to those obtained bv lUmetsu et alj d201 ll) . who did not include the 
2-halo term in their analysis, and provides a consistency check of 
our numerical scheme. 

The third best fit is the Stadel profile with Xv = 0.522 and 
X = 0.223 ± 0.040, followed by the Einasto profile with Xv = 0.602 
and n = 4.31 ±0.75. Next, the Sersic profile resulted in Xv = 0.663 
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Figure 4. Best fit curves and %v f° r tne phenomenological models stud- 
ied, ne glecting the 2-halo term, along with data points from lUmetsu et all 
1201 ll) . The bottom panel shows the relative difference between the best fits 
of each model and NFW. 
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Table 3. Fit results for the phenomenological models neglecting the 2-halo 
term. Columns defined as in Table[T] 



with n = 2. 69 ±0.41. Finally, for the Hernquist profile, with 2 scale 
parameters like NFW, we obtained Xv = 0.706. These results are 
summarized in Table [3] 

Fig. [5] shows the fits for the theoretical models when we ne- 
glect the 2-halo term. The best fit comes from the DARKexp model, 
for which %y = 0.598. The best fit value for the shape parameter 
was 90 = 3. 24 ±0.48. 

Among the Isothermal Sphere and its variants, the best fit is 
for the former, which provides %y = 2. 195 with X = (5.62 ± 0. 15) x 
10 4 . For the "Kang & He" model, with y = 3/5, we find x v = 2.265 
with X = (5.477 ± 0.003) x 10 4 and = 6.1 ± 0.6. Finally, for the 
"Kang & He 2" model with y = 4/5, we obtain %l = 2.385 with 
X= (5.52 ± 0.36) x 10 4 and// = 75 ±237. The results of the fits of 
these theoretical models are summarized in Table [4] 

As can be seen in Figs.|4]and[5] almost all the best fit density 
profiles remain below the data points in the outer regions. The ex- 
ceptions are the Isothermal Sphere variants, whose fits nonetheless 
fail badly. This trend shows the need for including the 2-halo term 
in the analysis. 




Figure 5. Best fit curves and x v for the theoretic al models stud i ed, ne glect- 
ing the 2-halo term, along with data points from lUmetsu et all J201 lh . The 
bottom panel shows the relative difference between the best fits of each 
model and DARKexp model. 
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Table 4. Fit results for the theoretical models neglecting the 2-halo term. 
Columns defined as in TablepT] 

7 DISCUSSION 

We have used observed data for the surface mass density of 4 clus- 
ters of similar mass and redshift to study various models for cluster 
density profiles. Under the assumption that the stacked data pro- 
vide a fair representation of the mean radial density profile of dark 
matter halos, we can include effects of large-scale structure at large 
radii and investigate how appropriately each model describes the 
average properties of clusters of this mass and redshift, and directly 
compare models against each other. 

For the phenomenological profiles, we have found that the 
Stadel model provides the best fit, followed by gNFW and Einasto. 
Nonetheless, the performance of all these profiles, including the 
standard NFW, is very similar, as is often the case in numerical sim - 
ulations l lGao et al.ll2008l ; lMerritt et al.l2005l : lNavarro et al.ll201fj) . 

For the theoretically motivated profiles, both the Isothermal 
Sphere and the Kang & He models give poor results compared to 
the phenomenological models. This can be attributed in part to the 
fact that both models produce cored density profiles, while the sim- 
ulations and the data used here favor cuspy profiles, or cored just in 
the innermost, inaccessible region, like for the Einasto and Stadel 
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Figure 6. Best fit curves fo r Stadel, DARKexp, Einasto and NFW models, 
along with data points from lUmetsu et all \20\ ll) . The bottom panel shows 
differences relative to the Stadel profile. 



profiles. Moreover, the outer region is not well described by these 
models, which behave like r~ 2 , while the data favor a behavior 
closer to r~ 3 . 

The best theoretical fit to data is obtained with the DARKexp 
model. This model provides an excellent fit, even compared with 
the performance of the phenomenological profiles. This is interest- 
ing, since this model has a dynamical basis justification. We did 
not inves tigate the role of the velocity anisotropy in this model. 
However, I Williams et all J2010h have compared the model to sim- 
ulated data and shown that the typical anisotropy profiles do not 
alter significantly the predicted density profile and that the DARK- 
exp model is a better match to the Einasto profile than to the NFW 
profile. This is in agreement with our findings as seen in Fig. [6] 
which shows our best fits for the Stadel, DARKexp, Einasto and 
NFW profiles. 

It is interesting to note that the Einasto profile with n ~ 6 
generally fits ACDM simulation s better tha n NFW, as has been 
noticed by iNavarro etall J2004I) . Moreover iMamon et"ail d2O10T) 
found n ~ 5 for dark matter halos of hydrodynamical cosmological 
simulations, while in this work we have found n ~ 3.8 . Final ly, by 
fitting rotation curves of spiral galaxies, IChemin et al.l d20 1 lb have 
found best fits for the Einasto profile with n as low or lower than 
unity. This sequence may suggest a correlation between the param- 
eter n and the importance of the gas component in the dynamics of 
the system. 

Concluding, we find that, while the observed data is best fit by 
phenomenological models, there is a similarly good fit by the theo- 
retical DARKexp model. We find that including the 2-halo term in 
the analysis is important and the best-fit ranking somewhat change 
if it is neglected. Should the agreement with observational data hold 
for clusters observed at different ranges of mass and redshift, one 
could argue that theoretical models such as DARKexp may provide 
a dynamical basis for the observed dark matter density profiles. 
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